library(vegan)
speabu<-read.csv('abund.csv',header=TRUE, row.names=1)
#transform data (spe abundance is high skewed and large values that are valid but highly influenctial) - doing a log transformation but add 1 since log zero is undefined
speabu.log<-log(speabu[,-1]+1)
#select a matrix coefficient - here doing a Bray-Curtis
speabu.d<-vegdist(speabu.log,"bray")
#perform PCoA. k is the number of primical compoents that should be extracted from this distance matrix (max here is min(col,row)-1), eig return eigenvalues default is false, if add is TRUE a constant is added to each value in dissim matrix so that resulting eigenvalues are non-negative (default here is FALSE)
spe.pcoa<-cmdscale(speabu.d,k=35,eig=TRUE,add=TRUE)
#look at it
spe.pcoa
#want to know the % variation explained by any principal component? then divide each eigenvalue by the sum of the eigenvalues across all PCs. Below looking at first 5 PCs
spe.pcoa$eig[1:5]/sum(spe.pcoa$eig)*100
#look for signifiance of each PC by compare eigenvalues to expectations according to broken stick model
plot(spe.pcoa$eig/sum(spe.pcoa$eig)*100-bstick(45)*100,xlab="PC", ylab="Actual-random percent variation explained") 
abline(h=0)
#look at ordination plot of PC and PC2
ordiplot(spe.pcoa, choices=c(1,2), type="text", display="sites", xlab="PC1 (15%)", ylab="PC2 (11%)")
#or maybe want to look at PC 1 and 3
ordiplot(spe.pcoa, choices=c(1,3), type="text", display="sites", xlab="PC1 (15%)", ylab="PC3 (9%)")
#if interested in looking at 1st 3 PCs simultaneously
pairs(spe.pcoa$points[,1:3],pch=24,col="red")
#Now calculate the PCA loadings (i.e. variable weights)
#here performs a linear regression between each of the original desriptors (species)and the scores from each PC. A permutation test is used to assess statistical significance
vec.sp<-envfit(spe.pcoa$points, speabu.log, perm=1000)
#another way to do the same thing is to extract first 2 PCs from PCoA object using the scores() function
vec.sp<-envfit(scores(spe.pcoa), speabu.log, perm=1000)
#look at values
vec.sp
#look closer at whats in vec.sp
names(vec.sp)
names(vec.sp$vectors)
#now plot the eigenvectors on the ordi plot
ordiplot(spe.pcoa, choices=c(1,2), type="text", display="sites", xlab='PC1 (15%)', ylab="PC2 (11%)")
plot(vec.sp, p.max=0.01, col="blue")